      program AMRID
C=====================================================================================
c
c     Original written by: Rosalinda de Fainchtein
c
c     Edited by: Yao, Cheng @ Harbin Institute of Techonology
c     2012-3-28
c
c     Reference:
c     Lohner, R."An Adaptive Finite Element Scheme for Transient Problems in CFD",
c     Comp. Meth. Appl. Mech. Eng.,1987,61, pp.323-338.
c
c     Available in NASA:
C     http://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19960016746_1996037138.pdf
C=====================================================================================
      parameter(mpoin=40000,melem=80000,period=0.D0)
      integer intmat(melem,2)
      real xp(mpoin),ce(melem,2),rmatm(mpoin),rlen(melem)
      real rho(mpoin),rhov(mpoin),rhoE(mpoin)
      real p(mpoin),v(mpoin),flux(mpoin,3)
      real diff(melem,3)
      real vmax(mpoin)
      real du(melem,3), dup(mpoin,3)

      integer mrhist(melem,5),lrefe(melem),ldere(melem)
      integer ipact(mpoin)
      integer iptemp(mpoin),ietemp(melem)
      integer l(mpoin)
      real tempea(melem),tempec(melem),errore(melem)
      real temppa(mpoin),temppb(mpoin),temppc(mpoin),errorp(mpoin)

c
c-----Get initial conditions and parameters
c

      call init(gamma,cour,ntime,di,ctore,ctode,epsil,nrmax,
     &          nbuff,ntref,npoin,nelem,mpoin,melem,period,
     &          intmat, mrhist, lrefe, ldere, ipact,
     &          xp, rho, p, v, rhov, rhoE,
     &          ibdry_l, ibdry_r )

      gamma1=gamma-1.D0
c
c-----Calculate element lengths, inverse mass matrix, and
c     coefficients for evaluating x derivatives on the elements.

      call geom(mpoin, melem, npoin, nelem, intmat, xp, ce, rmatm,
     &          rlen, period, mrhist, ipact )

c
c-----Refine a few times before starting
c

      do iloop=1,5

c
c-----Refinement Routines
c
c-----Computethe error function at each element and mark elements
c that require refinement or de-refinement

        call grerror(nelem, npoin, melem, mpoin, rho,
     &               tempea, tempec, temppa, temppb, temppc,
     &               epsil, ctore, ctode, nbuff,errorp,
     &               errore, mrhist, ipact,
     &               lrefe, ldere, ce, rlen, intmat,
     &               ietemp,iptemp,ibdry_l,ibdry_r )

c
c-----Coarsen elements marked for derefinement
c


        call unref(melem, mpoin, nelem, npoin, intmat,
     &             mrhist, ldere,ipact )

c
c-----Refine elements marked for refinement
c


        call refine(melem, mpoin, nelem, npoin, intmat, nrmax,
     &              mrhist, xp, rho, rhov, rhoE, lrefe, ipact)

c
c-----Recompute element lengths, inverse mass matrix and coefficients
c     for computing x derivatives within each element

        call geom(mpoin, melem, npoin, nelem, intmat, xp, ce, rmatm,
     &            rlen, period, mrhist, ipact )
c
c-----End of Refinement Routines
c

      end do
c
c-----BEGINTHEMAINLOOPOVERTIME STEPS-----
c==============================
      do i=1,ntime

c
c-----Refinement Routines
c

        if (mod(i,ntref) .eq. 0 ) then
c       
c-----  Compute the error function at each element and mark elements
c       that require refinement or de-refinement
c

          call grerror(nelem, npoin, melem, mpoin, rho,
     &                 tempea, tempec, temppa, temppb, temppc,
     &                 epsil, ctore, ctode, nbuff, errorp,
     &                 errore, mrhist, ipact,
     &                 lrefe, ldere, ce, rlen, intmat,
     &                 ietemp,iptemp,ibdry_l,ibdry_r )
        
c       
c-----   Coarsen elements marked for derefinement
c       
        
          call unref(melem, mpoin, nelem, npoin, intmat,
     &               mrhist, ldere,ipact )
c       
c-----   Refine elements marked for refinement
c       
        
          call refine(melem, mpoin, nelem, npoin, intmat, nrmax,
     &                mrhist, xp, rho, rhov, rhoE, lrefe, ipact)
c       
c-----  Re-compute element lengths, inverse massmatrix and coefficients
c       for computing x derivatives within each element
        
          call geom(mpoin, melem, npoin, nelem, intmat, xp, ce, rmatm,
     &              rlen, period, mrhist, ipact )
        
        end if
c
c-----End of Refinement Routines
c

c
c-----Computefluxes at the points and allowed time step.
c

        call fluxld(mpoin, melem, npoin, nelem, gamma,cour,
     &              intmat, rho, rhov, rhoE, p, v,
     &              rlen, vmax, flux, dt, mrhist, ipact )
        
c       
c-----  Computechanges in conserved variables on elements, and
c       scatter-add them to the point array dup
        
        call deltau(mpoin, melem, npoin, nelem,intmat,
     &              du,dup,flux,dt,ce,rlen, cour,di,
     &              rho,rhov,rhoE,diff,vmax,mrhist)
        
c       
c-----  Hard-Wall boundary conditions
c       
        
        if (period .eq. 0. ) then
          dup(ibdry_l,2)=0.
          dup(ibdry_r,2)=0.
        end if
        
c       
c-----  Advance the solution.
c       
        do ip=1,npoin
          if (ipact(ip) .eq. 1 ) then
          rho(ip)=rho(ip)-dup(ip,1)*rmatm(ip)*dt
          rhov(ip)=rhov(ip)-dup(ip,2)*rmatm(ip)*dt
          rhoE(ip)=rhoE(ip)-dup(ip,3)*rmatm(ip)*dt
          p(ip)= gamma1*(rhoE(ip)-0.5D0*rhov(ip)*v(ip))
          end if
        end do
c
c----- Output.
c      =======
c
        if ( mod(i, 10) .eq. 0 ) then
          write(50,*)' it',i,'	,dt=',dt
          write(50,*)' it npoin(active) nelem(active)'

c
c-----How many active points (neac)?
c

          npac=0
          do ip=1,mpoin
            if (ipact(ip) .eq. 1 ) then
            npac=npac+1
            end if
          end do

c
c-----How many active elements (neac)?
c

          neac=0
          do ie=1,melem
            if (mrhist(ie,5) .eq. 1 ) then
              neac=neac+1
              l(intmat(ie,1))=mrhist(ie,4)
            end if
          end do
          
          write(50,*) i,npac ,neac, npoin, nelem
          write(50,*)'     xp     rho     v',
     &     '     rhoE   p       I '
          
          do ip=1,npoin
            if (ipact(ip) .eq. 1 ) then
              write(50,6) xp(ip),rho(ip),v(ip),rhoE(ip),p(ip),l(ip)
            end if
          end do
        end if
      end do
c
c---END OF MAIN TIME STEP LOOP--
C
5     format(8(2x,f6.3))
6     format(5(2x,f6.3),1x,i1)
      stop
      end program


C*****************************************************************************
      subroutine init(gamma,cour,ntime,di,ctore,ctode,epsil,nrmax,
     &          nbuff,ntref,npoin, nelem, mpoin, melem,
     &          period, intmat, mrhist, lrefe, ldere, ipact,
     &          xp, rho, p, v, rhov, rhoE,
     &          ibdry_l, ibdry_r )
      real xp(mpoin),rho(mpoin),p(mpoin),v(mpoin),rhov(mpoin)
      real rhoE(mpoin)
      integer intmat(melem,2),mrhist(melem,5),lrefe(melem)
      integer ldere(melem),ipact(mpoin)
      integer ibdry_l,ibdry_r
      parameter(pi=3.141592654)
c
c-----Set constants, etc.
c

      gamma= 1.4D0
      gamma1=gamma-1.D0
      cour =0.5D0       ! Courant number
c      ntime=3000        ! Number of time steps
      ntime=300        ! Number of time steps
      di =0.5D0         ! Diffusion coefficient
      ctore=0.25D0      ! Threshold of error fcn. to refine
      ctode=0.1D0       ! Threshold to de-refine
      epsil=0.005D0     ! Refinement filter coefficient
      nbuff=3           ! No. of buffer layers
      ntref=2           ! No. of time steps per refinement
      nrmax=3           ! Max. No. of Refinement Levels.

C
C-----Set the initial conditions
C

c      npoin=128
      npoin=70
      nelem=npoin-1
      do ie=1,nelem
        intmat(ie,1)=ie
        intmat(ie,2)=ie+1
      end do
c
c-----boundary conditions.
c
      if (period .ne. 0.) then
        intmat(nelem,1)=nelem
        intmat(nelem,2)=1
      else
        ibdry_l=1
        ibdry_r=npoin
      end if


      do ie=1,melem
        do ir=1,5
          mrhist(ie,ir)=0
        end do
        lrefe(ie)=0
        ldere(ie)=0
      end do
      do ie=1,nelem
        mrhist(ie,5)=1
      end do


      do ip=1,npoin
        ipact(ip)=1
      end do

      do ip=npoin+1,mpoin
        ipact(ip)=0
      end do
C
C-----Initial Configuration for Sod Shock Tube Problem
C

      do ip=1,npoin/2
        rho(ip)=8.D0
        p(ip)=10.D0
      end do
      do ip=npoin/2,npoin
        rho(ip)=1.D0
        p(ip)=1.D0
      end do

      do ip=1,npoin
        xp(ip) =float(ip-1)*6.4D0/float(npoin)  !"float" is used to forcely change integer number to float number.
        v(ip) =0.d0
        rhov(ip)=rho(ip)*v(ip)
        rhoE(ip)=p(ip)/gamma1+0.5D0*rho(ip)*v(ip)*v(ip)
      end do

C
C-----Write out the initial conditions
C

      it=0
      write(50,*)'initial conditions'
      write(50,*)'      it      npoin      nelem'
      write(50,*) it,npoin,nelem
      write(50,*)'      xp      rho      v ',
     &     '      rhoE    p '
      do i=1,npoin
        write(50,6) xp(i),rho(i),v(i),rhoE(i),p(i)
      end do
6     format (5 (2x, e12.5) )
      return
      end subroutine


C*****************************************************************************
      subroutine geom(mpoin, melem, npoin, nelem, intmat, xp,
     &                ce, rmatm, rlen, period, mrhist, ipact )
C
C----- COMPUTE THE GEOMETRY ARRAYS FOR AID FE FLOW SIMULATION.
C
      integer mpoin,melem,npoin,nelem,intmat(melem,2)
      real xp(mpoin),ce(melem,2),rmatm(mpoin),rlen(melem)
      integer mrhist(melem,5),ipact(mpoin)

C
C-----Initialize the mass matrix
C

      do ip=1,npoin
        rmatm(ip)=0.D0
      end do

C
C-----Initialize element arrays
C

      do ie=1,nelem
      rlen(ie)=0.D0
      end do
C
C-----Compute the length of the elements (rlen)
C
      do ie=1,nelem
        if (mrhist(ie,5) .eq. 1 ) then
          rlen(ie)=xp(intmat(ie,2))-xp(intmat(ie,1))
        end if
      end do

c
c-----Periodic b.c.?
c
      if ( period .ne. 0.D0 ) then
        do ie=1,nelem
          if (abs(rlen(ie)) .gt. period/2. ) THEN
            rlen(ie)=period+rlen(ie)
          END IF
        end do
      end if

C
C-----The basis functions gradient
C

      do ie=1,nelem
        if (mrhist(ie,5) .eq. i ) then
          ce(ie,2)=1.D0/rlen(ie)
          ce(ie,1)=-1.D0*ce(ie,2)
        end if
      end do
c
c-----The inverse lumped mass matrix.
c

      do ie=1,nelem
        if (mrhist(ie,5) .eq. 1 ) then
          do in=1,2
            rmatm(intmat(ie,in))=rmatm(intmat(ie,in))+rlen(ie)
          end do
        end if
      end do
c
c-----Invert the sum
c

      do ip=1,npoin
        if ( ipact(ip) .eq. 1 ) then
          rmatm(ip)=2.D0/rmatm(ip)
        end if
      end do

      return
      end subroutine


C*****************************************************************************
      subroutine grerror(nelem, npoin, melem, mpoin, vref,
     &                   tempea, tempec, temppa, temppb, temppc,
     &                   epsil, ctore, ctode, nbuff, errorp,
     &                   errore, mrhist, ipact,
     &                   lrefe, ldere, ce, rlen, intmat,
     &                   ietemp,iptemp,ibdry_l,ibdry_r )
c
c-----THIS SUBROUTINE COMPUTES THE H-2 SEMINORM ERROR ESTIMATOR
c     [ errorp=temppa/ ( temppb+epsil*temppc) ]
c     (Ref: Lohner)

      integer intmat(melem,2),mrhist(melem,5),ipact(mpoin)
      integer lrefe(melem),ldere(melem)
      integer iptemp(mpoin),ietemp(melem)
      real ce(melem,2),rlen(melem)
      real tempea(melem),tempec(melem),errore(melem)
      real temppa(mpoin),temppb(mpoin),temppc(mpoin),errorp(mpoin)
      real vref(mpoin)
      integer ibdry_l,ibdry_r

c
c  nbuff=3
c
c----Initialize element arrays
c
      do ie=1,melem
        tempea(ie)=0.D0
        tempec(ie)=0.D0
        errore(ie)=0.D0
        lrefe(ie) =0.D0
        ldere(ie) =0.D0
      end do
C
C-----Initialize point arrays
C

      do ip=1,mpoin
        temppa(ip)=0.D0
        temppb(ip)=0.D0
        temppc(ip)=0.D0
        errorp(ip)=0.D0
      end do
C
C-----Sums at the elements.
C
      do ie=1,nelem
        if (mrhist(ie,5) .eq. 1) then
          do in=1,2
            tempea(ie)=tempea(ie)+ce(ie,in)*vref(intmat(ie,in))
            tempec(ie)=tempec(ie)+abs(ce(ie,in) )*
     &             abs(vref(intmat(ie,in)) )
          end do
        end if
      end do
c
c-----Scatter to the points.
c
      do ie=1,nelem
        if (mrhist(ie,5) .eq. 1) then
          do in=1,2
            temppa(intmat(ie,in))=temppa(intmat(ie,in)) +
     &                            rlen(ie)*ce(ie,in)*tempea(ie)
            temppb(intmat(ie,in))=temppb(intmat(ie,in)) +
     &                            rlen(ie)*abs(ce(ie,in)) *
     &                            abs(tempea(ie))
            temppc(intmat(ie,in))=temppc(intmat(ie,in)) +
     &                            rlen(ie)*abs(ce(ie,in)) *
     &                            tempec(ie)
          end do
        end if
      end do
C
C-----Boundary terms
C
      temppa(ibdry_l)=0.D0
      temppa(ibdry_r)=0.D0
C
C-----The error estimate at the points
C

      do ip=1,npoin
        if (ipact(ip) .eq. 1 ) THEN
          errorp(ip)=abs(temppa(ip))/(temppb(ip)+epsil*temppc(ip) )
        END IF
      end do

C
C-----Compute the error estimate at the elements (max. of node values).
C
      do ie=1,nelem
        if (mrhist(ie,5) .eq. 1) THEN
          errore(ie)=amax1(errorp(intmat(ie,1)),errorp(intmat(ie,2)))
        END IF
      end do
C
C-----Mark the elements for refinement or de-refinement.
C
      do ie=1,nelem
        if (mrhist(ie,5).Eq. 1) then
          if (errore(ie).Gt. ctore) lrefe(ie)=1
          if (errore(ie).Lt. ctode) ldere(ie)=1
        end if
      end do
c
c-----Buffer Layers
c
      do ib=1,nbuff

        do ip=1,npoin
          iptemp(ip)=0
        end do

        do ie=1,nelem
          if (mrhist(ie,5).eq. 1) then
            do in=1,2
              if (lrefe(ie) .eq. 1) then
                iptemp(intmat(ie,in))=1
              end if
            end do
          end if
        end do
      
        do ie=1,nelem
          if (mrhist(ie,5).eq. 1) then
            do in=1,2
              if (iptemp(intmat(ie,in)) .eq. 1) then
                lrefe(ie)=1
              end if
            end do
          end if
        end do

      end do
c
c-----Space Continuity on levels of refinement (only necessary on 1D)
c
      do iloop=1,5
      
        do ie=1,nelem
          ietemp(ie)=0
          if (mrhist(ie,5).eq.l) then
            ietemp(ie)=mrhist(ie,4)+lrefe(ie)
          end if
        end do
      
        do ip=1,npoin
          iptemp(ip)=0 
        end do
      
        do ie=1,nelem
          if (mrhist(ie,5).eq. 1) then
            do in=1,2
              if (ietemp(ie) .gt. iptemp(intmat(ie,in)) )then
                iptemp(intmat(ie,in))=ietemp(ie)
              end if
            end do
          end if
        end do
      
        do ie=1,nelem
          if (mrhist(ie,5).eq. 1) then
            do in=1,2
              if (iptemp(intmat(ie,in))-mrhist(ie,4) .gt. 1) then
                lrefe(ie)=1
              end if
              if (iptemp(intmat(ie,in))-mrhist(ie,4) .ge. i) then
                ldere(ie)=0
              end if
            end do

            ietemp(ie)=mrhist(ie,4)+lrefe(ie)
            do in=1,2
              if (ietemp(ie) .gt. iptemp(intmat(ie,in)) ) then
                iptemp(intmat(ie,in))=ietemp(ie)
              end if
            end do
          end if
        end do !End of ie

      end do   !End of iloop
      return
      end subroutine 


C*****************************************************************************
      subroutine unref(melem, mpoin, nelem, npoin, intmat,
     &                 mrhist, ldere,ipact )
c
c-----THIS ROUTINE DE-REFINES THOSE ELEMENTS MARKED F0R DE-REFINEMENT
c
      integer intmat(melem,2)
      integer mrhist(melem,5),ldere(melem)
      integer ipact(mpoin)

c-----AMR arrays
c     ==========

c     mrhist(ie,1) == parent element number. (default=0)
c               2) == first child element number. (default=0)
c               3) == second child element number. (default=0)
c               4) == level of refinement. (default=0)
c               5) == element active? ( 1=yes )
c     ldere(ie) =1 =======> de-refine
c               =0 =======> do not de-refine
c
c-----Loop over the elements
c
      do ie=1,nelem
        if ( ldere(ie).eq. 1 .and. mrhist(ie,4).ne. 0 ) then
c
c-----Identify elements involved
c
        iparent=mrhist(ie,1)
        ichildl=mrhist(iparent,2)
        ichild2=mrhist(iparent,3)
c
c-----Derefine iff both children are marked for de-refinement (4 steps)
c
          if ( ldere(ichildl).eq. 1 .and. ldere(ichild2).eq. 1 ) then
c           
c-----(1) Identify and _"deactivate" points not assoc, w/parent.
c           
            
            do in=1,2
              ipl=intmat(ichildl,in)
              ip2=intmat(ichild2,in)
              if(ipl .ne. intmat(iparent,1) .and. 
     &           ipl .ne. intmat(iparent,2)) then
                ipact(ipl)=0
              end if
            
              if(ip2 .ne. intmat(iparent,1) .and.
     &           ip2 .ne. intmat(iparent,2)) then
                ipact(ip2)=0
              end if
            end do
            
c           
c-----(2) Disable children elements
c           
            mrhist(ichildl,5)=0
            mrhist(ichild2,5)=0
c           
c-----(3) Enable parent element 
c           
            mrhist(iparent,5)=1
c           
c-----(4) Unmarkthe children elements
c           
            ldere(ichildl)=0
            ldere(ichild2)=0
          end if  ! (   if all children marked)
        end if    ! (   if one child marked)
      end do      ! (end loop over elements)
      
      return
      end subroutine


C*****************************************************************************
      subroutine refine(melem, mpoin, nelem, npoin, intmat, nrmax,
     &                  mrhist, xp, rho, rhov, rhoE, lrefe, ipact )

c
c-----THIS ROUTINE REFINES TH0SE ELEMENTS MARKED F0K REFINEMENT
c
      integer intmat(melem,2)
      real xp(mpoin)
      real rho(mpoin),rhov(mpoin),rhoE(mpoin)
      integer mrhist(melem,5),lrefe(melem)
      integer ipact(mpoin)
C
C-----AMR arrays
C     ==========
C     mrhist(ie,1) == parent element number. (default=0)
C               2) == first child element number. (default=0)
C               3) == second child element number. (default=0)
C               4) == level of refinement. (default=0)
C               5) == element active? ( l=yes )
C     lrefe(ie) =1 ======> refine
C               =0 ======> do not refine
C
C----- Loop over the elements
C

      nrmaxl=nrmax-1
      do ie=1,nelem
        if (lrefe(ie).eq. 1 .and. mrhist(ie,4).le. nrmaxl ) then
          npoin=npoin+1 !new point
          ipact(npoin)=1
          if (npoin .gt. mpoin) then
            write(*,*)'Please increase mpoin. Needed:',npoin
            stop
          end if
          nelem2=nelem+2
          if (nelem2 .gt. melem) then
            write(*,*)'Please increase melem. Needed:',nelem2
            stop
          end if
      
          xp(npoin)  =0.5D0*(xp(intmat(ie,1))+xp(intmat(ie,2)))
          rho(npoin) =0.5D0*(rho(intmat(ie,1))+rho(intmat(ie,2)))
          rhov(npoin)=0.5D0*(rhov(intmat(ie,1))+rhov(intmat(ie,2)))
          rhoE(npoin)=0.5D0*(rhoE(intmat(ie,1))+rhoE(intmat(ie,2)))

C
C-----New Element Connectivity (Input new elements at the end of the list)
C
          intmat(nelem+1,1)=intmat(ie,1)
          intmat(nelem+1,2)=npoin
          intmat(nelem+2,1)=npoin
          intmat(nelem+2,2)=intmat(ie,2)
C
C-----Disable the parent elements and enable the children
C
          mrhist(ie,5)=0
          mrhist(nelem+1,5)=1
          mrhist(nelem+2,5)=1
C
C-----Modify the rest of the mesh refinement history parameters
C
          mrhist(ie,2)=nelem+1
          mrhist(ie,3)=nelem+2
          mrhist(nelem+1,1)=ie
          mrhist(nelem+2,1)=ie
          mrhist(nelem+1,4)=mrhist(ie,4)+1
          mrhist(nelem+2,4)=mrhist(ie,4)+1
C
C-----Modify the current number of elements
C
          nelem=nelem+2
        end if
      end do
      return
      end


C*****************************************************************************
      subroutine fluxld(mpoin, melem, npoin, nelem, gamma, cour,
     &                  intmat, rho, rhov, rhoE, p, v,
     &                  rlen, vmax, flux, dt, mrhist, ipact )
C
C-----Computethe fluxes at the points and the time step.
C     ==================================================

      integer mpoin,melem,npoin,nelem,intmat(melem,2)
      real rho(mpoin),rhov(mpoin),rhoE(mpoin)
      real p(mpoin),v(mpoin),flux(mpoin,3)
      real rlen(melem),vmax(mpoin)
      integer mrhist(melem,5),ipact(mpoin)
      cbig=1.E30
      gamma1=gamma-1.d0

      do ip=1 ,npoin
        if ( ipact(ip).eq. 1 ) then
          v(ip) =rhov(ip)/rho(ip)
          p(ip)= gamma1 * (rhoE(ip)-0.5D0*rhov(ip)*v(ip) )
          flux(ip,1)=rhov(ip)
          flux(ip,2)=rhov(ip)*v(ip)+p(ip)
          flux(ip,3)=rhoE(ip)*v(ip)+p(ip)*v(ip)
        end if
      end do
C
C-----TIME STEP: dt.
C     =============

c
c-----Compute the "Cmaximum velocity" (sound speed+v) at the points.
c

      do ip=1,npoin
        if ( ipact(ip).eq. 1 ) then
          vmax(ip)=sqrt(gamma*p(ip)/rho(ip))+abs(v(ip))
        end if
      end do
c
c-----and the time step...
c
      dt= cbig
      do ie=1,nelem
        if (mrhist(ie,5).eq. 1 ) then
          do in=1,2
            dt=amin1(dt,rlen(ie)/vmax(intmat(ie,in)))
          end do
        end if
      end do

      dt=cour*dt

      return
      end subroutine


C*****************************************************************************
      subroutine deltau(mpoin, melem, npoin,nelem,intmat,
     &                  du,dup,flux,dt,ce,rlen,cour,di,
     &                  rho,rhov,rhoE,diff,vmax,mrhist)
C
c-----Compute the du's and dup's
c

      integer mpoin,melem,npoin,nelem, intmat(melem,2)
      real du(melem,3), dup(mpoin,3),flux(mpoin,3),dt
      real ce(melem,2),rlen(melem)
      real rho(mpoin),rhov(mpoin),rhoE(mpoin),diff(melem,3)
      real vmax(mpoin)
      integer mrhist(melem,5)
      cbig=1.e20
      csmall=-1.e20
c
c-----Initialize du
c

      do im=1,3
        do ie=l,nelem
          du(ie,im) =0.D0
          diff(ie,im)=0.D0
        end do
      end do
C
C-----The minimumdelta x.
C

      rlmin=cbig
      rlmax=csmall
      do ie=1,nelem
        if (mrhist(ie,5).eq. 1 .and. rlen(ie).Lt. rlmin) THEN
          rlmin=rlen(ie)
        END IF

        if (mrhist(ie,5).eq. 1 .and. rlen(ie).gt. rlmax) THEN
          rlmax=rlen(ie)
        END IF
      end do
C
C-----Initialize dup
C

      do im=1,3
        do ip=1,npoin
          dup(ip,im)=0.D0
        end do
      end do
C
C-----BEGIN LOOP OVER THE ELEMENTS
C

      do ie=1,nelem
        if (mrhist(ie,5) .eq. 1 ) then
C       
C-----Compute du at the elements
C         
          do im=1,3
            do in=1,2
              du(ie,im)=du(ie,im)+ce(ie,in)*flux(intmat(ie,in),im)
            end do
          end do
C         
C-----Computethe diffusion terms at the elements
C         

          vv=amax1(vmax(intmat(ie,1)),vmax(intmat(ie,2)))
          do in=1,2
            diff(ie,1)=diff(ie,1)+ce(ie,in)*rho(intmat(ie,in))
     &                  *vv
            diff(ie,2)=diff(ie,2)+ce(ie,in)*rhov(intmat(ie,in))
     &                  *vv
            diff(ie,3)=diff(ie,3)+ce(ie,in)*rhoE(intmat(ie,in))
     &                  *vv
          end do
          dfact=0.5D0*di
          do im=1,3
            diff(ie,im)= dfact *diff (ie,im)*rlen(ie) *rlmin
          end do
c         
c-----Scatter to the points.
c         
          do im=1,3
            do in=1,2
              dup(intmat(ie,in),im)=dup(intmat(ie,in),im)    +
     &                               0.5D0*rlen(ie)*du(ie,im)+
     &                               ce(ie,in)*diff(ie,im)
            end do
          end do
        end if
      END DO

C
C-----END OF THE ELEMENT LOOP--
c
      return
      end subroutine